{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Velocity Data Download Script for Validation\n",
    "#### Working Script to pull large velocity datasets and later compare to ICESAT2-derived velocities\n",
    "\n",
    "ICESat-2 hackweek  \n",
    "June 15, 2020  \n",
    "Lynn Kaluzienski"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os,re,h5py\n",
    "import requests\n",
    "import zipfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import pyproj\n",
    "import scipy, sys, os, pyproj, glob\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import Point, Polygon\n",
    "import pointCollection as pc\n",
    "\n",
    "# Import some of the scripts that we have written\n",
    "import sys\n",
    "sys.path.append(\"/home/jovyan/surface_velocity/scripts\")\n",
    "from loading_scripts import atl06_to_dict\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-86 -81 -81 -86 -86]\n",
      "[-55 -55 -65 -65 -55]\n"
     ]
    }
   ],
   "source": [
    "#From Ben Smith's code loading in .tif file, running into issues likely with directories\n",
    "data_root='/srv/shared/surface_velocity/FIS_Velocity/'\n",
    "#spatial_extent = np.array([-102, -76, -98, -74.5])\n",
    "spatial_extent = np.array([-65, -86, -55, -81])\n",
    "lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "print(lat)\n",
    "print(lon)\n",
    "# project the coordinates to Antarctic polar stereographic\n",
    "xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# get the bounds of the projected coordinates \n",
    "XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "#Originally tried to load data from a local directory, should change to shared directory\n",
    "Measures_vx=pc.grid.data().from_geotif(os.path.join(data_root,'Measures2_FIS_Vx.tif'), bounds=[XR, YR])\n",
    "Measures_vy=pc.grid.data().from_geotif(os.path.join(data_root,'Measures2_FIS_Vy.tif'), bounds=[XR, YR])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9fa19425fadf4f2d82238cadbce4f787",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'viridis', 'clim': [-100, 100], 'extent': array([-887878.65809562, -400567.81840096,  200333.90920048,\n",
      "        561654.87344315]), 'origin': 'lower'}\n",
      "{'cmap': 'viridis', 'clim': [-100, 100], 'extent': array([-887878.65809562, -400567.81840096,  200333.90920048,\n",
      "        561654.87344315]), 'origin': 'lower'}\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8bdcb64bed3547069d9b9dbb5ad097b1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load a line and plot\n",
    "\n",
    "data_root='/srv/shared/surface_velocity/'\n",
    "field_dict={None:['delta_time','latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "\n",
    "rgt = \"0848\"\n",
    "cycle=\"03\"\n",
    "filename = glob.glob(os.path.join(data_root, 'FIS_ATL06', f'*ATL06_*_{rgt}{cycle}*_003*.h5'))[0]\n",
    "\n",
    "D=atl06_to_dict(filename,'/gt2l', field_dict=field_dict, index=None, epsg=3031)\n",
    "\n",
    "# show the velocity map:\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(121)\n",
    "Measures_vx.show(cmap='viridis', clim=[-100,100])\n",
    "plt.plot(xy[0,:], xy[1,:],'k')\n",
    "plt.title('Measures X-Velocity')\n",
    "plt.plot(D['x'],D['y'],'r')\n",
    "plt.gca().set_aspect('equal')\n",
    "\n",
    "plt.subplot(122)\n",
    "Measures_vy.show(cmap='viridis', clim=[-100,100])\n",
    "plt.plot(xy[0,:], xy[1,:],'k')\n",
    "plt.title('Measures Y-Velocity')\n",
    "plt.plot(D['x'],D['y'],'r')\n",
    "plt.gca().set_aspect('equal')\n",
    "\n",
    "#plt.tight_layout()\n",
    "\n",
    "# Interpolate the Measures velocities along the line and plot (note that these are speeds for now)\n",
    "\n",
    "vx = Measures_vx.interp(D['x'],D['y'])\n",
    "vy = Measures_vy.interp(D['x'],D['y'])\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(121)\n",
    "plt.plot(D['x_atc'],vx)\n",
    "plt.title('X-Velocity')\n",
    "plt.subplot(122)\n",
    "plt.plot(D['x_atc'],vy)\n",
    "plt.title('Y-Velocity')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# add measures along track to dictionary\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lenght of the first segment is  19.47695252813686  this is close to the teoretical elngth of 20 m\n",
      "The angle of the track for vel transfromation is:  0.8401541611474018  radians\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# should be close to 20\n",
    "hypoth = np.sqrt((D['x'][0]-D['x'][1])**2+(D['y'][0]-D['y'][1])**2)\n",
    "print(\"Lenght of the first segment is \", hypoth , \" this is close to the teoretical elngth of 20 m\")\n",
    "\n",
    "adjacent = abs(D['x'][0]-D['x'][1])\n",
    "\n",
    "angle = np.arccos(adjacent/hypoth)\n",
    "\n",
    "print(\"The angle of the track for vel transfromation is: \" ,angle, \" radians\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "vx = Measures_vx.interp(D['x'],D['y'])\n",
    "vy = Measures_vy.interp(D['x'],D['y'])\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
